function K = bending_bb_2(dofs,V,T,d,Int_UU_cell,Int_UV_cell,Int_UW_cell,Int_VV_cell,Int_VW_cell,Int_WW_cell);
%        K = bending_bb_2(dofs,V,T,nt,d,areas);
% This function computes the bending matrix $\int_\Omega \Delta \phi_i \Delta \phi_j$. 
% This matlab program is copyrighted @2001 by Ming-Jun Lai and Paul Wenston
% through University of Georgia Research Foundation, Inc..
% %comput the total striff matrix
if nargin == 4    
% prepare the global bending matrix for saving time
    d_max = max(d);
    Int_UU_cell = cell(d_max,1);
    Int_UV_cell = cell(d_max,1);
    Int_UW_cell = cell(d_max,1);
    Int_VV_cell = cell(d_max,1);
    Int_VW_cell = cell(d_max,1);
    Int_WW_cell = cell(d_max,1);
    for degree = 2:d_max
        Mat = build_dim2(degree-2,degree-2);
        Id = diag(ones((degree+1)*(degree+2)/2,1));
        Du = dirder(Id,1,0,-1); % the direction derivate on v1-v3
        Duu = dirder(Du,1,0,-1); % the direction derivate on v1-v3
        Duv = dirder(Du,0,1,-1);  % blending partial derivatives
        Dv = dirder(Id,0,1,-1);  % the direction derivate on v2-v3
        Dvv = dirder(Dv,0,1,-1); % the direction derivate on v2-v3
        Int_UU_cell{degree} = Duu'*Mat*Duu;
        Int_UV_cell{degree} = Duu'*Mat*Dvv;
        Int_UW_cell{degree} = Duu'*Mat*Duv;
        Int_VV_cell{degree} = Dvv'*Mat*Dvv;
        Int_VW_cell{degree} = Dvv'*Mat*Duv;
        Int_WW_cell{degree} = Duv'*Mat*Duv;
    end
end

nt = size(T,1);

% get the areas of all triangles.
x21 = V(T(:,2),1) - V(T(:,1),1);
x31 = V(T(:,3),1) - V(T(:,1),1);
y21 = V(T(:,2),2) - V(T(:,1),2);
y31 = V(T(:,3),2) - V(T(:,1),2);
areas = (x21.*y31 - x31.*y21)/2;

m = 0;
for k=1:nt
    md = (d(k)+1)*(d(k)+2)/2;
    m = m + md^2;
end
Indx1 = zeros(m,1);
Indx2 = zeros(m,1);
S = zeros(m,1);
pos_start = 1;
for k = 1:nt
    Int_UU = Int_UU_cell{d(k)};
    Int_UV = Int_UV_cell{d(k)};
    Int_UW = Int_UW_cell{d(k)};
    Int_VV = Int_VV_cell{d(k)};
    Int_VW = Int_VW_cell{d(k)};
    Int_WW = Int_WW_cell{d(k)};
    V1=V(T(k,1),:);V2=V(T(k,2),:);V3=V(T(k,3),:);
    u = (V2(2)-V3(2))^2+(V2(1)-V3(1))^2;
    v = (V1(2)-V3(2))^2+(V1(1)-V3(1))^2;
    w = 2*((V2(2)-V3(2))*(V1(2)-V3(2))+(V2(1)-V3(1))*(V1(1)-V3(1)));
    LocK = (u*u*Int_UU - u*w*(Int_UW+Int_UW') + v*v*Int_VV +...
        w*w*Int_WW - v*w*(Int_VW+Int_VW') + u*v*(Int_UV+Int_UV'))/(16*areas(k)^3);   % J = 2*tri_area(K)
    [i,j,s] = find(LocK);
    loc_dof = dofs(k,1):dofs(k,2);
    L = length(i);
    Indx1(pos_start:(pos_start + L-1)) = loc_dof(i);
    Indx2(pos_start:(pos_start + L-1)) = loc_dof(j);
    S(pos_start:(pos_start + L-1)) = s;
    pos_start = pos_start + L;   
end
dim = max(max(dofs));
K = sparse(Indx1(1:(pos_start-1)),Indx2(1:(pos_start-1)),S(1:(pos_start-1)),dim,dim);